Implementation of the Pn-approximation for Radiative Heat Transfer on Openfoam

This work presents an OpenFOAM implementation of the PN approximation for radiative heat transfer, including higher orders P3, P5, and P7. Also described is a procedure which enables the sequential numerical computations of the coupled partial differential equations (PDEs) by re-expressing the boundary conditions in matrix form so that individual boundary conditions can be associated with each PDE. The implementation of the software programs are verified with derived analytical solutions for 1-D slabs with constant and variable properties, and are also tested with various orientations in order to demonstrate the geometric invariance properties of the 3-dimensional PN formulation. A few examples taken from the literature are also considered in this work and could be taken as benchmark solutions for the PN approximations. INTRODUCTION The radiative transfer equation (RTE) is an integrodifferential equation in seven independent variables (3 space and 2 directional, temperature, and wavenumber) [1], which is exceedingly difficult to solve. As a result, approximate solution methods to the RTEs such as the spherical harmonics method (SHM), discrete ordinates method (DOM), the finite volume method (FVM), or Monte Carlo method are frequently employed to solve radiation problems. Each of these approximating methods have their well-known advantages and disadvantages. The DOM/FVM method has the advantage that it is easy to implement and can use the finite volume discretizations which are used ∗ASME Life Fellow; Address all correspondence to this author. Phone: (209) 228-4113. in computational fluid dynamics (CFD) computations. Higher accuracy can be obtained using the DOM/FVM by increasing angular discretizations. The disadvantages of DOM/FVM is that results can become inaccurate when there are large aspect ratios, geometries with sharp corners, or if the extinction coefficient varies by orders of magnitude in the computational domain [2,3]. The PMC method is an exact method with statistical error that decreases with larger number of photon bundles. The PMC is often used to benchmark other RTE solvers when exact solutions are unavailable. However, PMC can be much slower than other RTE solvers. The SHM, which is the focus of this paper, offers an approximate solution to the radiative transfer equation (RTE) by transforming the RTE to a system of elliptic PDEs. This method approximates the radiative intensity as a truncated series of spherical harmonics that decouple the directional and spatial variation of the intensity field. For optically thick cases, the lowest order of the SHM, P1, can provide acceptable accuracy for the intensity field. However, for optically thin to intermediate cases, the accuracy of P1 is known to diminish [1]. Despite the shortcomings of P1, applications of the higher order SHM have been rather limited, due to exponentially increasing mathematical complexity. Modest and Yang [4] formulated a methodology to reduce the (N +1)2 first-order PDEs to a set of N(N +1)/2 second-order, elliptic PDEs. Their work also included the formulation of the N(N +1)/2 Marshak boundary conditions. The methodology was designed to be applicable for arbitrary three dimensional geometries. Further development of the PN-approximation for isotropic scattering and its Marshak boundary conditions was recently presented in [5]. The newer formulation for isotropic scattering, reduces the complexity of the SHM formulation because many of the coefficient tensors involved can be combined to obtain all1 Copyright c © 2013 by ASME symmetric operators in the set of PN equations. There are only a few literature examples with applications of high order SHM. Yang and Modest [6] demonstrated P3 on a two dimensional rectangular enclosure with an isotropically scattering medium with a hot wall segment, and also a two-dimensional triangular enclosure with an absorbing-emitting medium with variable absorption coefficient. Ravishankar et al. [3] demonstrated P3 in a steady state homogenous combustion of a methane-air mixture under laminar flow conditions for optical thicknesses τL = 0.25 and 0.5. Modest [5] applied P30 and P5 to a square enclosure with an absorbing-emitting medium with strongly variable temperature and absorption coefficients. The absorption coefficients were adjusted in order to analyze optically thin, intermediate, and thick situations. The work presented here employs the full three-dimensional PN approximation formulated in [5], where the only simplification is to neglect anisotropic scattering. The set of governing differential equations consist of N(N + 1)/2 unknowns coupled by a set of elliptical PDEs that are complex because of several second-order mixed derivatives that are involved. These mixed derivatives makes it challenging to obtain analytical solutions for simple results in 2and 3-D. However, for simple 1-D cases, analytical solutions are readily obtained. The 1-D problems are important as they serve to verify that the PN implementation accurately employs each of the terms that are involved. Additionally, one can verify that the resulting incident radiation profiles computed by PN are invariant with respect to coordinate rotations. The work presented here also tests 2and 3-D examples. The 2-D example includes the sample problem provided in [5], which is further extended to demonstrate P7. A version of this 2-D example is also applied to a 3-D cylinder case. With these example problems, the purpose of this work is to 1. Outline useful procedures, such as manipulation of the boundary conditions, in order to implement the SHM in the OpenFOAM CFD software package [7], 2. verify consistency of solutions of PN solutions for geometrically invariant problems, 3. demonstrate the relative accuracy of higher order PNapproximations on 1-, 2-, and 3-D examples, 4. and illustrate PN solutions for axisymmetric cases. The procedures introduced here take advantage of the SHM formulation for isotropic scattering [5]. Similar procedures are more involved for the anisotropic scattering formulation and are not discussed here. PN FORMULATION Governing Differential Equations The PN approximation is based on representing the intensity field I(τ,s) as a series products of intensity coefficients functions and of spherical harmonics Ym n , whereby the spatial (~τ = ∫ βd~x, β is the extinction coefficient) and the directional (ŝ) dependencies are decoupled,